Simulating the dynamics of dispersal and dispersal ability in fragmented populations with mate‐finding Allee effects

Abstract We consider the spatial propagation and genetic evolution of model populations comprising multiple subpopulations, each distinguished by its own characteristic dispersal rate. Mate finding is modeled in accord with the assumption that reproduction is based on random encounters between pairs of individuals, so that the frequency of interbreeding between two subpopulations is proportional to the product of local population densities of each. The resulting nonlinear growth term produces an Allee effect, whereby reproduction rates are lower in sparsely populated areas; the distribution of dispersal rates that evolves is then highly dependent upon the population's initial spatial distribution. In a series of numerical test cases, we consider how these dynamics affect lattice‐like arrangements of population fragments, and investigate how a population's initial fragmentation determines the dispersal rates that evolve as a habitat is colonized. First, we consider a case where initial population fragments coincide with habitat islands, within which death rates differ from those that apply outside; the presence of inhospitable exterior regions exaggerates Allee effect‐driven reductions in dispersal ability. We then examine how greater distances separating adjacent population fragments lead to more severe reductions in dispersal ability. For populations of a fixed initial magnitude, fragmentation into smaller, denser patches leads not only to greater losses of dispersal ability, but also helps ensure the population's long‐term persistence, emphasizing the trade‐offs between the benefits and risks of rapid dispersal under Allee effects. Next, simulations of well‐established populations disrupted by localized depopulation events illustrate how mate‐finding Allee effects and spatial heterogeneity can drive a population's dispersal ability to evolve either downward or upward depending on conditions, highlighting a qualitative distinction between population fragmentation and habitat heterogeneity. A final test case compares populations that are fragmented across multiple scales, demonstrating how differences in the relative scales of micro‐ and macro‐level fragmentation can lead to qualitatively different evolutionary outcomes.

based on random encounters between pairs of individuals, so that the frequency of interbreeding between two subpopulations is proportional to the product of local population densities of each. The resulting nonlinear growth term produces an Allee effect, whereby reproduction rates are lower in sparsely populated areas; the distribution of dispersal rates that evolves is then highly dependent upon the population's initial spatial distribution. In a series of numerical test cases, we consider how these dynamics affect lattice-like arrangements of population fragments, and investigate how a population's initial fragmentation determines the dispersal rates that evolve as a habitat is colonized. First, we consider a case where initial population fragments coincide with habitat islands, within which death rates differ from those that apply outside; the presence of inhospitable exterior regions exaggerates Allee effect-driven reductions in dispersal ability. We then examine how greater distances separating adjacent population fragments lead to more severe reductions in dispersal ability. For populations of a fixed initial magnitude, fragmentation into smaller, denser patches leads not only to greater losses of dispersal ability, but also helps ensure the population's long-term persistence, emphasizing the trade-offs between the benefits and risks of rapid dispersal under Allee effects. Next, simulations of well-established populations disrupted by localized depopulation events illustrate how mate-finding Allee effects and spatial heterogeneity can drive a population's dispersal ability to evolve either downward or upward depending on conditions, highlighting a qualitative distinction between population fragmentation and habitat heterogeneity. A final test case compares populations that are fragmented across multiple scales, demonstrating how differences in the relative scales of micro-and macro-level fragmentation can lead to qualitatively different evolutionary outcomes.

K E Y W O R D S
Allee effect, dispersal, evolution of dispersal ability, fragmented population, mate finding, reaction-diffusion model 1 | INTRODUC TI ON

| Fragmentation and dispersal ability
Changes to a population's habitat can lead to changes in the phenotypic traits that the population exhibits. This is perhaps most visible through phenomena such as island gigantism and island dwarfism (Benítez-López et al., 2021;Lomolino, 2005;McClain et al., 2013;Raia & Meiri, 2006;Van Valen, 1973). Aside from changes in size, insular populations can also evolve drastic morphological differences, such as those that lead to flightlessness in birds (Wright et al., 2016) or increased woodiness in plants (Lens et al., 2013). It has been suggested that most (if not all) of these evolutionary changes reflect the same consistent trend toward reduced dispersal ability on islands, or even beyond islands (Filin & Ziv, 2004;Lomolino, 2005;Waters et al., 2020;Whittaker & Fernández-Palacios, 2007). Darwin famously pondered insular dispersal ability loss by considering an analogy with shipwrecked mariners facing a choice between clinging to the shipwreck or swimming away (Lomolino, 2009). Individuals with an innate tendency to stay put (i.e., slower dispersers) would remain, while those predisposed to swim away (i.e., faster dispersers) would leave. In this way, slower dispersers could come to consolidate themselves there, reducing the local population's dispersal ability. A variety of mathematical analyses repeatedly predicted similar tendencies toward dispersal ability loss, even beyond island habitats (Asmussen, 1983;Balkau & Feldman, 1973;Filin & Ziv, 2004;Hastings, 1983;Holt, 1985;Johnson & Gaines, 1990). The downward evolution of dispersal ability in these models is typically driven by the adverse effects suffered disproportionately by rapid dispersers as they traverse harmful features of their habitats, such as dangerous boundaries or gradients in environmental quality. In this way, these models typically assume that some form of environmental heterogeneity, rather than spatial isolation, is the primary factor shaping the evolution of dispersal characteristics. In doing so, they often draw conclusions which do not depend on the initial distribution of a population throughout its environment.
Elsewhere, spatial isolation has often been used to explain these insular phenomena, whether they occur on true islands bounded by water, or on habitat islands, where other forms of isolation lead to "island effects" such as body size change and dispersal ability loss (Amburgey et al., 2021;Cayuela et al., 2019;Haila, 2002;Incagnone et al., 2015;Lens et al., 2013;McClain et al., 2006;Merckx et al., 2018). The isolation of population fragments from one another has long been recognized to have a complex influence on the evolution of populations and species beyond islands as well (Kisel & Barraclough, 2010;Losos et al., 2010;MacArthur & Wilson, 2001;Whittaker et al., 2017;Whittaker & Fernández-Palacios, 2007). Spatial isolation can alter the selective pressures that shape the evolution of a segment of a population (Jessop et al., 2018;Lomolino et al., 2012;McClain et al., 2013;Millien, 2004) by sheltering it from predators, competitors, and environmental hazards, or conversely by separating individuals from resources or from their own peers, including potential mates. Even within Darwin's "shipwreck" scenario, the prospects of survival for slower or faster dispersers would depend on the details of the shipwreck, including how the mariners and ship fragments were distributed through space relative to one another and to various features of the surrounding environment. In realistic populations and habitats too, spatial heterogeneity can result in a complex interplay between environments, populations, and the patterns of dispersal and dispersal ability that evolve.
Both habitats and the populations that inhabit them can simultaneously exhibit heterogeneity, patchiness, or fragmentation, with the distances between population fragments recognized as a factor in the likelihood of successful dispersal (Bowler & Benton, 2005;Conradt et al., 2000). Empirical observations of genetic rescue-by which migration-driven gene flow reintroduces genetic diversity into isolated population fragments, supporting their continued adaptation and survival-demonstrate how flows between separate fragments can be mediated by individuals with high dispersal ability (Bell et al., 2019;Ingvarsson, 2001;Räsänen & Hendry, 2008;Whiteley et al., 2015). Natural or anthropogenic disturbances can alter habitats while also affecting population densities across affected areas, fragmenting habitats and populations; variations in body size and dispersal ability have been observed to follow these events (Brisson et al., 2003;Griffiths & Brook, 2014;Merckx et al., 2018;Palkovacs et al., 2012). When populations expand into unpopulated areas, dispersal abilities have been observed to evolve upward along the advancing edges of population fragments (Bénichou et al., 2012;Bouin et al., 2012;Deforet et al., 2019;Holt et al., 2004;Hughes et al., 2007;Léotard et al., 2009;Phillips et al., 2010;Travis et al., 2009). These examples highlight the potentially crucial role of heterogeneity of a population's distribution throughout an environment-population fragmentation, as distinct from habitat fragmentation-in shaping the evolution of dispersalrelated traits. This study thus focuses upon the lesser studied aspects of population fragmentation, in terms of which issues of spatial isolation can be disentangled from heterogeneity in the underlying habitat. Using a reaction-diffusion model, we demonstrate how the details of a population's initial fragmentation in space can play an important role in shaping the dispersal characteristics that develop when a genetically diverse population reproduces sexually. Before presenting the model, we will review related previous work on reaction-diffusion and metapopulation models of coupled spatialgenetic dispersal dynamics, as well as mate-finding Allee effects, by which spatial aggregation, rather than isolation, becomes advantageous for sexually reproducing populations.

T A X O N O M Y C L A S S I F I C A T I O N
Biogeography, Evolutionary ecology, Landscape ecology, Movement ecology, Population ecology, Population genetics, Quantitative genetics, Spatial ecology 1.2 | Modeling dispersal and the evolution of dispersal ability

| Reaction-diffusion models
A population's spatial movements can change its patterns of exposure to its environment, while also affecting how frequently different types of individuals within the population interact with one another. These changes, in turn, alter the birth and death processes that shape the population's genetics, including those traits that determine how it moves through space. This results in a complex feedback between dispersal and reproduction that can be readily captured by reaction-diffusion equations. Skellam's (1951) pioneering reaction-diffusion model describes how a population, represented by a density function, evolves under simultaneous processes of random-walk dispersal and reproduction (see Box 1, Equation B1).
Considering a model habitat encircled by a "zone of absolute extinction," Skellam concluded that a population with a slower dispersal rate would grow more quickly, while a population of faster dispersers would grow more slowly, or even decay, as it spilled more rapidly outward into the habitat's lethal exterior. Filin and Ziv (2004) later invoked this result to explain the apparent universal tendency toward dispersal ability loss on islands: subpopulations with slower dispersal rates would grow faster than subpopulations of rapid dispersers. However, this explanation relies on the assumption of passive dispersal across a lethal "absorbing" island boundary. Its heuristic arguments also overlook the possibility that subpopulations distinguished by different dispersal rates can mate, interacting through reproduction to potentially "rescue" one another from extinction. These limitations demonstrate the need for reactiondiffusion models that can (1) accommodate a wider range of domains and boundaries, and (2) more explicitly account for interactions between subpopulations with different dispersal rates.
A number of studies sought to further develop the pioneering work of Skellam and others (e.g., Kierstead & Slobodkin, 1953) by applying reaction-diffusion models to investigate ecological problems in greater depth (Britton, 1986;Ōkubo et al., 2001); a comprehensive review is given by Cantrell and Cosner (2004). These analyses considered reaction-diffusion dynamics on domains with more general shapes and boundary conditions, while sometimes also accommodating spatial heterogeneity among per-capita growth rates or dispersal rates within a domain's interior (see Box 1, Equations B2 and B3). These analyses provided a more thorough theoretical understanding of how a habitat's shape, boundaries, and interior source-sink dynamics can affect a dispersing population's long-term persistence, making predictions about the critical patch sizes required for survival (Cantrell & Cosner, 2001, 2004. These insights were used to formulate more generalized reaction-diffusion approaches toward island biogeography (Cantrell et al., 1996;Cantrell & Cosner, 1994, 2001, while remaining applicable to a wider variety of scenarios of interest in landscape ecology. Meanwhile, other reaction-diffusion modeling efforts explicitly modeled the interactions between coexisting subpopulations with different dispersal traits. Unlike models that focused on the longterm persistence of populations of individuals all sharing the same dispersal rate, Dockery et al. (1998) explicitly modeled variability among dispersal abilities. Their approach considered the coevolution of multiple population density functions, each representing a phenotype characterized by its own distinct dispersal rate (see Box 1, Equation B4), and coupled to the other phenotypes through competition for resources and small mutations. By tracking how the relative abundances of slower and faster dispersers would change as they dispersed through an environment with spatially varying carrying capacity, the model predicted a universal tendency toward dispersal ability loss (Dockery et al., 1998). Other approaches have since obtained similar results using models formulated with continuous rather than discretized dispersal rates (Lam & Lou, 2017), and ongoing research has continued to use reaction-diffusion approaches that incorporate environmental heterogeneity in new ways (Cantrell et al., 2020;Wickman et al., 2017).

| Allee effects
An Allee effect (Courchamp et al., 2008) operates when, at lower population densities, increasing density has a positive effect on fitness and reproduction rates. Aggregation, rather than isolation, becomes advantageous. A strong Allee effect applies when the effect can go beyond merely slowing growth rates to cause a net population decline. A variety of mechanisms produce Allee effects; for example, spatial aggregation by animals can facilitate cooperation in hunting, foraging, or defense, while in plants, higher vegetation density can help maintain favorable soil conditions that support further growth (Rietkerk et al., 2004). Mate-finding Allee effects specifically associated with sexual reproduction can arise when, for example, animals in sparsely populated areas seldom encounter potential mates, or as pollen propagated by plants into sparsely populated areas too often fails to reach conspecifics (Davis et al., 2004).
In the context of reaction-diffusion models, per-capita growth rates (see Box 1) can depend on local population density, with the appropriate mathematical form of density dependence determined by the specific mechanisms at hand (Aronson & Weinberger, 1978;Cantrell et al., 1996;Du et al., 2019;Du & Shi, 2007;Shi & Shivaji, 2006). Allee effects can arise from mate-finding processes BOX 1 Related reaction-diffusion models. In Skellam's (1951) seminal reaction-diffusion model, the evolution of a population density function Ψ(x, y) is described by a partial differential equation that combines simultaneous processes of random-walk dispersal (with characteristic step size proportional to a constant a) and reproduction (at a constant per-capita growth rate c > 0): y 2 is the Laplacian operator. Skellam considered Equation (B1) on a circular domain of radius r 0 , beyond the outer edge of which lie a "zone of absolute extinction," while enforcing an "absorbing" (i.e., Dirichlet) boundary condition (Pudjaprasetya, 2018) of Ψ = 0 along the edge for continuity. Regardless of the population's initial distribution throughout the domain, solutions have a dominant mode consisting of a dome-shaped density function that grows (or decays) exponentially at a spatially uniform rate k = c − a 2 j 2 1 ∕ 4r 0 (where j 1 ≈ 2.405). Populations with faster dispersal rates would traverse the "absorbing" boundary in greater numbers, inhibiting their growth; this lethal boundary effect becomes more exaggerated on smaller domains. For a given domain size r 0 , populations with excessive dispersal rates (a > 2 √ cr 0 ∕ j 1 ) fail to maintain densities sufficient to support net growth, and so decay to extinction; similarly, for a population with a given dispersal rate a, there is a critical patch size (r 0 > a 2 j 2 1 ∕ (4c)) below which the population cannot persist. More general reaction-diffusion models have taken forms such as where the functions a(x, y) and c(x, y) can accommodate spatial variations among dispersal rates and per-capita growth rates, respectively (Cantrell & Cosner, 2004). Boundary conditions were also formulated much more generally, allowing for a hybrid of (partial) absorption and (partial) reflection: where � ⃗ n is an outward normal vector with respect to the domain boundary, and the function (x, y) describes the fraction of individuals impinging upon the boundary at (x, y) that can traverse it (Cantrell & Cosner, 2004). Treatment of eigenvalue problems that arise from Equations (B2) and (B3) yielded more general conclusions about the overall rates of population growth or decay: the rate of population loss across a boundary is proportional to the largest eigenvalue of the diffusion operator on the domain Ω, which is determined primarily by the patch area and not by irregularities in the shape of the boundary Ω. Even when "no-flux" (i.e., Neumann) boundary conditions (Pudjaprasetya, 2018) are applied, analyses of source-sink dynamics within the domain lead to predictions regarding critical patch sizes and long-term persistence (Cantrell & Cosner, 2001, 2004. The tendency toward reduced dispersal ability in spatially-heterogeneous habitats was studied within a reaction-diffusion framework by Dockery et al. (1998), who considered the coevolution of multiple population density functions Ψ k (for k = 1, … , n), each representing a distinct phenotype k with characteristic dispersal rate a k , in environments with spatially varying carrying capacities K(x, y): where > 0 is small and M ki encodes the relative frequencies of random mutations from phenotype i into phenotype k. Spatial heterogeneity in the environment (i.e., in K(x, y)) was shown to shift the relative abundances of the phenotypes in favor of slower dispersers.
In contrast with this study, in which subpopulations interact directly via sexual reproduction, the functions Ψ k in that case were coupled to one another only via the logistic growth term, which depends on total population density. The small linear mutation term was used to confirm the robustness of the model's conclusions with respect to mutations.
due to the relative rarity of encounters between potential mates in sparsely populated areas (Boukal & Berec, 2002Gascoigne et al., 2009;Lutscher et al., 2022;McCarthy, 1997); these often share a mathematical form anticipated by Volterra and later termed a bimolecular collision model (Dennis, 1989). In these models. the frequencies of encounters between members of different subpopulations are assumed to be proportional to the product of their densities. Some models separately account for male and female subpopulations by allowing for fluctuating sex ratios (Boukal & Berec, 2002Gascoigne et al., 2009 2 | ME THODS

| Dynamical equations
We consider the coevolution of n population density functions k (X, Y), each representing a distinct genotype k distinguished by its own characteristic dispersal ability a k (related to the distance traveled per unit time in random-walk movements). These dispersal abilities take one of n evenly spaced values, a k = a 0 k, for k = 1, … , n (where a 0 > 0). The environment is assumed to have a finite carrying capacity K, such that a logistic growth factor attenuates birth rates as the environment becomes saturated. The mate-finding process affects birth rates in accord with the assumption that births result from random encounters between pairs of individuals (as in bimolecular collision models; Dennis, 1989). By further assuming a constant, spatially uniform sex ratio, the local probability of encounters between members of subpopulations i and j becomes proportional to the product i j . A factor ij k describes the probability ( ∑ n k=1 ij k = 1 ) that parents of genotypes i and j will produce offspring of type k .
Deaths are assumed to occur randomly with probability d within each time increment. The population density function k representing genotype k thus evolves with respect to time according to the dynamical equation Y 2 is the Laplacian operator and all parameters are positive.
We focus on the case where offspring have an equal 50% chance of inheriting the genotype k of either parent, so that is then identical across all genotypes k, depending only upon the overall population density . As in some previous models (e.g., Dockery et al., 1998), then, any changes in the relative abundances of the different genotypes can be attributed unambiguously to their different dispersal rates. In contrast to those models, the mate-finding process modeled here results in a different nonlinear dependence of the per-capita growth rate upon total population .
Specifically, Equation (4) recalls the class of strong Allee effect growth terms studied previously for populations with a uniform dispersal rate (e.g., Amarasekare, 1998;Cantrell et al., 1996;Dennis, 1989;Du & Shi, 2007;Wang et al., 2011Wang et al., , 2019. These growth terms exhibit bistable "explosion/extinction" behavior (Du & Shi, 2007;Shi & Shivaji, 2006;Wang et al., 2011), always evolving toward one of the two possible outcomes: (1) successful colonization of the environ- Beyond addressing questions of long-term population persistence, though, the inclusion of multiple dispersal genotypes here enables us to consider the distributions of dispersal ability that evolve from different fragmented initial conditions.
If population densities are expressed as fractions of carrying where ∇ 2 = 2 x 2 + 2 y 2 now, and Ψ = ∑ n i=1 Ψ i is the rescaled overall population density. Two dimensionless parameters remain: a rescaled birth parameter ≡ K r 0 ∕a 0 2 b and rescaled death rate ≡ r 0 ∕a 0 2 d. This illustrates that for given a set of initial population configurations, the task of exploring the range of dynamics possible

| Numerical scheme
We simulate the dynamics of Equation (5) using a finite difference method (Pudjaprasetya, 2018), approximating the Laplacian operator ∇ 2 using a 9-point stencil (LeVeque, 2007) while applying periodic boundary conditions. We consider a computational grid with node spacing Δx spanning horizontal coordinates X = x 1 , … , x N x and vertical coordinates Y = y 1 , … , y N y . Discretized population density Simulation results are then summarized in terms of the total , and the population's mean dispersal ability, a(t) = 1 Unless otherwise noted, simulations were terminated when the value (t) ∕ Ψ 1 (t) (where is the standard deviation of genotype k = 1 density values Ψ 1 x i , y j ; t , and Ψ 1 is their mean) receded to below 10 −3 (i.e., when the slowest class of dispersers have nearly achieved a spatially uniform population density), or when the total overall population P(t) receded to below 10 −4 (extinction). Simulation code is available at https://osf.io/qy5gf/ ? view_only=9d069 efcd7 6e437 9a8a6 874b2 7dd2e4d.

| Fragmented population configurations
Rigorous studies of reaction-diffusion models have explored how patch/fragment geometry determines long-term outcomes by delineating the ranges of patch sizes, densities, or spacings over which survival or extinction will result (Cantrell & Cosner, 2001, 2004. While a similar analytical approach is beyond our scope, this study also systematically explores how different aspects of a population's spatial configuration-fragment sizes, densities, and spacingsaffect evolutionary outcomes under mate-finding Allee effects. To this end, we deal with idealized fragmented populations of which the characteristic sizes, densities, and spacings of fragments can be varied ( Figure 1). We detail the layouts of these configurations in the following.
To approximate an equilateral triangular lattice of circular patches with nearest-neighbor spacing s, grid points x i , y j are de- For each initial configuration, simulations can be repeated over a range of birth parameter and death rate values to explore how outcomes are affected by environmental conditions. This study comprises five test cases (summarized in Table 1). The first three test cases model sparse, fragmented populations dispersing within an otherwise unpopulated domain ( ext = 0). The extent of fragmentation is varied from trial to trial (in terms of spacing s in Experiments A and B, and in terms of radius r in Experiment C). The mean dispersal rates a that evolve are then observed for those populations that persist. In Experiment D, the layout is inverted to simulate an otherwise-saturated domain ( ext = K) in which circular regions are initially unpopulated ( int = 0); this can be seen as representing a well-established population following its disturbance by some spatially localized, catastrophic depopulation events. Experiment E considers populations that are fragmented across multiple spatial scales, with fragments forming a roughly self-similar lattice of lattices ( Figure 1c). Additional details about each of these test cases are discussed alongside simulation results below.

| Experiment A: Dispersal of population fragments between habitat islands
Numerous studies have identified gene flow between population fragments, including cases of "genetic rescue" between habitat islands, as important factors in evolution. Situations like these can be modeled as an "archipelago" of habitat patches separated from one another by regions with inhospitable-but neither strictly lethal nor impenetrable-conditions. Other work has applied reactiondiffusion approaches to investigate related issues of island biogeography or other complex habitats, but has not typically focused on how the initial arrangements of population fragments-as distinct from habitat fragments-might affect the extent of the subsequent changes in dispersal ability. Experiment A uses this "habitat islands" scenario to clarify and distinguish the potential roles of habitat and population fragmentation in complex scenarios like these. Its results provide context for the spatially homogeneous test cases that follow.
Circular population fragments, initially populated with density int (Figure 1a), are set to coincide with circular habitat islands wherein (x, y) = int , with initially unpopulated exterior regions where (x, y) = ext . While holding patch radii r, initial interior and exterior densities int and ext , and parameters and int constant across all trials (see Table 1), we repeat simulations over a range of values of patch-exterior death rate ext . The final mean dispersal rates a shown in Figure 2a are the values achieved when the condition 1 is first satisfied for all k; these do not represent steady states, but rather the states attained soon after the habitat has become saturated and Allee effects have ceased to play a primary role.
The curve representing the homogeneous environment case Over the longer run, however, spatial gradients in the death rate (x, y) will continue to drive a net flux of dispersers into the more lethal regions. Genotypes that disperse more quickly into these less hospitable regions will be disproportionately affected, gradually draining the population of its more rapid dispersers. The mean dispersal abilities that initially result from Allee effects (Figure 2a) will not persist in the long run. When habitat and population fragmentation coincide, the heterogeneity-driven dispersal ability loss observed in so many previous models will indeed occur here. However, since these habitat-driven changes can be orders of magnitude slower than those that result from Allee effects, the dispersal traits that initially evolve due to fragmentation can endure for a relatively

| Experiment B: The role of separation between population fragments
Experiment B retains the initial population configurations used in Experiment A, but now considers their dispersal within spatially homogeneous environments over a range of combinations of birth and death parameters. While holding the initial sizes, shapes, and densities of these patches consistent across all trials, we run a series of simulations that vary the spacing s separating the fragment centers. Lower values of s give more closely packed fragments; as s is increased, these fragments become increasingly isolated from one another. As simulations proceed from these fragmented initial conditions (t 1 in Figure 3), faster dispersers propagate first into these unpopulated regions, and so are more severely affected by Allee effects; meanwhile, slower dispersers remain densely aggregated around fragment centers, where they maintain higher growth rates (t 2 in Figure 3). When fragments are more closely spaced, dispersers from adjacent fragments more quickly meet and establish higher population densities in the gaps between patches (t 3 in Figure 3).
This reduces the losses of rapid dispersers suffered early on, retaining more of their genes within a population that (given sufficiently hospitable conditions) then proceeds to saturate the environment (t 4 in Figure 3). However, when fragments are separated by larger spacings s, the greater times taken to traverse the intervening gaps results in greater losses of rapid dispersers before overall population numbers stabilize. Once the population has established itself in greater numbers, surviving rapid dispersers may regain some advantage as they are able to occupy the sparsely populated regions between fragments first (t 3 − t 4 in Figure 3c,f). By initially dominating these regions, rapid dispersers can partially recover the dispersal ability lost during more precarious stages of evolution (t 2 − t 3 in Figure 3).
Even as a population saturates the environment, and its constituent genotypes become uniformly intermixed in space, the trauma of the initial fragmentation nonetheless remains "frozen in" to the population, having determined the relative abundances of different genotypes. The mean dispersal ability values attained decrease as the initial spacing s increases (Figure 2b). If the initial spacing s is further increased, then Allee effects become insurmountable and the entire population decays. Under more favorable environmental conditions (i.e., higher ratios ∕ ), Allee effect dispersal ability losses (and, eventually, extinctions) begin to take effect at larger spacings s ; the more favorable the conditions, the greater the amount of initial spatial isolation that can be overcome.
Animations for the simulations summarized in Figure 3 are shown in Movies S1 and S2, and animations for the simulations of Experiment B, representing the range of parameter values of s, , and included in Figure 2b, are accessible at https://osf.io/2qn8u/ ? view_only=7aa42 a040f 0b417 aa571 266a0 e465a9d.

| Experiment C: The role of fragment size and density
In Experiment B, the spacing s between fragment centers was used to track the effects of spatial separation upon the dispersal abilities that evolved. Increasing the patch spacing s in this way also results in an increased domain area (Figure 1a), which decreases the population's initial overall density (i.e., the ratio of total overall population to domain area (≈ P ∕ � √ 3s 2 � )). Experiment C aims to disentangle the geometric aspects of fragmentation from variations to the overall population density by varying patch size r in tandem with density int , so that the geometry of the initial configuration is altered without affecting its initial overall density.
Many results show significant reductions in mean dispersal ability a, with dynamics qualitatively similar to those observed in TA B L E 1 Initial conditions and parameter settings/ranges for the five test cases. Independent variables (i.e., quantities plotted along the horizontal axes in Figure 2) are italicized; Ranges in parentheses (e.g., "(350-400)") indicate the limits of the parameter regimes over which simulation trials were repeated. Entries for Experiments A through D correspond to the layout of Figure 1a. Entries for Experiment E refer to the (i) coarse fragmentation and (ii) fine fragmentation layouts displayed in Figure 1c, with S and R indicating the spacing and radius, respectively, of the "macro"-lattice of hexagonal fragments, with s and r describing the spacings and radii of the circular "micro"-fragments.

Experiment
Lattice spacing s Radius r Experiment B (Figure 2c). The observed decline in mean dispersal rates as initial patch radii r are decreased (i.e., progressing leftward

Birth parameter(s) Death rate(s)
in Figure 2c) corresponds to the decline observed for increasing spacing s in Experiment B: in both cases (as s is increased, or r is decreased), the unpopulated gaps between fragments are expanded.
Here, however, the non-monotonic dependence of the final mean dispersal rate a upon r (Figure 2c) illustrates the trade-off between high population density within fragments (which enhances a fragmented population's ability overcome Allee effects during the precarious initial stages of evolution) and smaller intervening distances dispersal ability (or extinction) as they struggle to populate the vast surrounding regions. On the other hand, when fragments are wider but less densely aggregated, populations may be spread too thin to overcome Allee effects even within fragments. Between these extremes, the advantages of closer spacings between fragments (which support the survival of faster dispersers) and higher densities within fragments (where slower dispersers stay put and reproduce) complement one another. All else being equal, an intermediate de- gree of spatial fragmentation here enhances the chances of survival compared to more-or less-fragmented initial configurations, including for more rapid dispersers who disproportionately suffer Allee effects at both extremes.

| Experiment E: Fragmentation across multiple spatial scales
The idealized fragmented configurations of the previous test cases demonstrate how the spatial separations, densities, and sizes of identical fragments affect the extent of dispersal ability losses (or gains) that evolve due to mate-finding Allee effects. More realistic populations can contain fragmentation at multiple spatial scales.
Experiment E considers more elaborate configurations where fragments are grouped into larger macro-fragments, which themselves are also arranged in a triangular lattice (Figure 1c). We consider two different population configurations that differ only in terms of their spatial fragmentation at the microscale. Both feature large hexagonal regions with "radii" of around 5 and lattice spacings of 13.2, and the same total initial population. In the first configuration, these macrofragments are more coarsely subdivided (with r = 1.0 and s = 2.2); in the second, they are more finely subdivided (with r = 0.5 and s = 1.1).
Simulating x,y,k Ψ k (x, y;t) ) at indicated times t 1 , t 2 , t 3 , and t 4 for spacing s = 2.15. (d) Time evolution of the total population of each genotype for spacing s = 2.40. (e) Time evolution of population mean dispersal ability a for spacing s = 2.40. (f) Relative population densities Ψ k (x, y;t) ∕ Ψ max (t) at indicated times t 1 , t 2 , t 3 , and t 4 for spacing s = 2.40. Note that by the times t 4 shown here, the relative abundances of each genotype have nearly reached their steady-state values; following t 4 , populations continue to diffuse, approaching spatially uniform distributions.

F I G U R E 4
Simulation results from Experiment E with = 385, = 37. (a) Time evolution of the total population of each genotype for the coarse fragmentation case. (b) Time evolution of population mean dispersal ability a for the coarse fragmentation case. (c) Relative population densities Ψ k (x, y;t) ∕ Ψ max (t) (where Ψ max (t) = max x,y,k Ψ k (x, y;t)) at indicated times t 1 , t 2 , t 3 , and t 4 for the coarse fragmentation case. (d) Time evolution of the total population of each genotype for the fine fragmentation case. (e) Time evolution of population mean dispersal ability a for the fine fragmentation case. (f) Relative population densities Ψ k (x, y;t) ∕ Ψ max (t) at indicated times t 1 , t 2 , t 3 , and t 4 for the fine fragmentation case. Note that by the times t 4 shown here, the relative abundances of each genotype have nearly reached their steady-state values; following t 4 , populations continue to diffuse, approaching spatially uniform distributions. regions, where rapid dispersers come to dominate slower dispersers (t 3 and t 4 in Figure 4a-c), leading to a net dispersal ability increase (Figure 4d-f). The effects of the population's initial microfragmentation are thus still felt even long after its traces are no longer visible in the population's spatial distribution (t 4 in Figure 4d).
These effects are manifest not just in the magnitude of dispersal change, but potentially also in its direction: upward or downward.
Full animations for the simulations summarized in Figure 4 are shown in Movies S3 and S4.

| DISCUSS ION
The model presented here combines several features of previous reaction-diffusion models in a novel way to study the evolution of dispersal by fragmented populations with an Allee effect. Its use of multiple interacting dispersal ability genotypes distinguishes it from other related reaction-diffusion models that included Allee effects (Cantrell et al., 1996;Shi & Shivaji, 2006). As in the work of Dockery et al. (1998), variability among dispersal abilities is modeled by discretizing the space of possible dispersal rates into several bins, each represented by its own characteristic dispersal ability (a k ), and then considering the coevolution of the population density functions (Ψ k ) that represent these subpopulations. In that model, direct interactions between different kinds of dispersers were not a prominent driver of evolution in their own right. Here, however, nonlinear interactions between these different genotypes play a more formative Reaction-diffusion models have long been used to demonstrate how-for populations confronted with heterogeneous habitats or hostile boundaries-slow dispersal can come to dominate over time.
The simulations presented here illustrate a case where multiple coexisting strategies-slower dispersers who remain anchored around dense population fragments, and rapid dispersers who colonize sparsely populated areas-can complement one another. Each contributes potential advantages at different phases of evolution; the relative magnitudes of these advantages depend on the distribution of the population in space. This sensitivity to initial conditions is another qualitative difference distinguishing the present model from a great deal of previous related work. Following Skellam (1951), many analytical reaction-diffusion approaches to dispersal have tended to focus only on the long-term outcomes approached after all traces of initial conditions have vanished. Our results, however, provide examples of highly path-dependent dispersal evolution, where transient dynamics determine the states that a population approaches in the long run. Even in homogeneous environments, where previously studied mechanisms eventually act to reduce dispersal ability, the consequences of these earlier stages of evolution can remain relevant over relatively long-time scales (Experiment A).
Within the current model, the details of this fragmentation can even determine the direction-downward or upward-in which dispersal ability evolves. This capacity for predicting upward evolution of dispersal rates under some conditions, even within an unchanging habitat, marks an interesting qualitative distinction of this model from the majority of related work. A few simple mechanisms here can lead either to downward or upward changes in dispersal ability depending on fragmentation patterns and environmental conditions. Given that the purportedly universal loss of dispersal ability on islands has been called into question even for the most passive of dispersers (Burns, 2018;García-Verdugo et al., 2017), models like this could be useful when navigating these issues. Some of the simulation results presented above exhibit features which recall various empirical results reported elsewhere, such as the upward evolution of dispersal ability along the expanding fronts of population fragments, which in some cases have been attributed to Allee effects (e.g., Tobin et al., 2007). Further work could explore the possibility of establishing more direct links between models and simulations like those of this study and more concrete empirical observations. These numerical experiments have explored the space of potential dynamics that can result from population fragmentation and Allee effects, demonstrating some qualitative differences of these dynamics from those of habitat fragmentation. A rigorous theoretical treatment of this model, more akin to that undertaken for other reactiondiffusion models (Cantrell & Cosner, 2004), could provide deeper insights. One advantage of this study's computational approach, though, is that simulation results allow for visualizations of the complex patterns of spatial propagation undergone by multiple interacting classes of dispersers alongside the population-level genetic changes in dispersal ability that result (as in Figures 3 and 4 and Movies S1-S4).
The real utility of such a model may lie in its application to more realistic and specific case studies, where the spatial patterns it predicts could be compared with observed spatial distributions. Previous work has applied related models to problems of conservation or pest control (e.g., Boukal & Berec, 2009;Du et al., 2019); in applications like these, simulations could be useful to help guide spatially targeted interventions. The modeling and simulation approach presented here, by providing visualizations of complex coupled processes of spatial propagation and genetic evolution, could help to facilitate these kinds of connections between models and applications in ways that more abstract, theoretical results alone might not.

ACK N OWLED G M ENTS
The authors wish to thank K. Christopher Beard and the University

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available from the Open Science Framework at https://osf.io/qy5gf/ ?view_ only=9d069 efcd7 6e437 9a8a6 874b2 7dd2e4d.